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Abstract. The main purpose of this article is to show how symmetry structures in par¬ 
tial differential equations can be preserved in a discrete world and reflected in difference 
schemes. Three different structure preserving discretizations of the Liouville equation are 
presented and then used to solve specific boundary value problems. The results are com¬ 
pared with exact solutions satisfying the same boundary conditions. All three discretizations 
are on four point lattices. One preserves linearizability of the equation, another the infinite¬ 
dimensional symmetry group as higher symmetries, the third one preserves the maximal 
finite-dimensional subgroup of the symmetry group as point symmetries. A 9-point in¬ 
variant scheme that gives a better approximation of the equation, but significantly worse 
numerical results for solutions is presented and discussed. 
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1 Introduction 

This article is part of a general program the aim of which is to make full use of the theory of 
Lie groups to study the solution space of discrete equations and in particular to solve difference 
equations [7, 8, 21, 24, 28, 29, 30, 42]. This is one of the areas to which Luc Vinet made 
important contributions [14, 15, 16, 27]. 

S. Lie introduced what is now called Lie groups as groups of transformations of the inde¬ 
pendent and dependent variables figuring in a system of differential equations [31, 36]. Of 
special importance are symmetry groups, transforming solutions into solutions. These may be 
point transformations, where new variables depend only on the old ones. They may be contact 
transformations, where the new variables depend also on the first derivatives of the dependent 
variables. They may also be generalized symmetries where the new variables can also depend 
on all derivatives of the old ones. 

*This paper is a contribution to the Special Issue on Exact Solvability and Symmetry Avatars in honour of 
Luc Vinet. The full collection is available at http://www.emis.de/journals/SIGMA/ESSA2014.html 
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Lie’s method is particularly powerful for ordinary differential equations (ODEs). A one¬ 
dimensional point symmetry group can be used to lower the order of the ODE by one. An n- 
dimensional (solvable) Lie point symmetry group can be used to decrease the order by n. Thus, 
if the ODE is of order n! < n Lie group theory can provide the general solution (i.e., one that 
satisfies arbitrary initial conditions) in explicit or implicit analytic form. Eor partial differential 
equations (PDEs) the Lie point symmetry group is used to decrease the number of independent 
variables in the equation and to provide special solutions (group invariant solutions), satisfying 
particularly symmetrical boundary conditions. 

The aim of this general program is to extend the use of Lie symmetry groups to difference 
systems (AS), i.e., to difference equations together with the lattice they are written on. 

The program has two complementary aspects, an analytical and a numerical one. 

The aim of the analytical aspect is to determine the maximal symmetry group of the AS, i.e., 
the group of transformations that takes solutions into solutions, and then to use it to obtain 
exact analytic solutions, at least special ones, if possible general ones. The AS to which the 
approach is applied can come from the study of discrete physical, chemical, biological or other 
systems, for which symmetries play an important role. Among them we mention phenomena in 
crystals, or in atomic or molecular chains. 

On the other hand AS can be obtained by discretizing ODEs, or PDEs, that have nontrivial 
symmetry groups reflecting fundamental physical laws such as Galilei, Lorentz, or conformal 
invariance. At the scale of the Planck length space-time may very well be discrete. In this case 
continuous equations are approximations (continuous limits) of discrete ones. Erom the physical 
point of view the symmetries are very important and should be preserved, e.g., when studying 
quantum field theories on lattices. 

One way of preserving symmetries in a discretization of continuous equations (the one used 
in this article) is to use symmetry adapted lattices that themselves transform under the group 
action. This greatly enlarges the set of equations for which symmetry preserving discretization 
is possible. We will however see that in some cases only a subgroup of the Lie point symmetry 
group can be preserved as point symmetries. 

The numerical aspect of our program is the following. When solving an ODE or PDE nu¬ 
merically it is always necessary to replace the continuous equation by a difference system. This 
can be done in a standard manner, applicable to all equations, simply by replacing derivatives 
by discrete derivatives. The other possibility takes us directly into the field of geometric integra¬ 
tion [20, 22, 33, 34]. The idea is to focus on some important feature of the underlying problem 
and to preserve it in the discretization. Such a feature may be, for instance linearizability, 
hamiltonian structure, integrability in the sense of the existence of a Lax pairs and generalized 
symmetries or point and contact symmetries. We are concentrating on point symmetries and 
exploring the possibility and usefulness of including them in numerical calculations. 

Earlier work has shown that for first-order ODEs preserving a one-dimensional symmetry 
group provides an exact discretization [40]. Eor second-order ODEs preserving a 3-dimensional 
symmetry group often provides analytically solvable schemes (either via a Lagrangian [11, 12] 
or via the adjoint equation method [9]). Eor third- and higher-order ODEs symmetry preserving 
discretization provides numerical solutions that are, usually, closer to exact ones then those 
obtained by other methods, specially near to the singularities [5, 39]. Eor previous work on 
PDEs see [2, 3, 4, 6, 10, 19, 25, 26, 37, 38, 41]. 

Several recent articles [1, 23, 38] were devoted to discretizations of the Liouville equation [32] 

^xy — 6 , (lA) 

or its algebraic version 
UUxy UxUy — U , 


u = e . 


( 1 . 2 ) 
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The Liouville equation is of interest for many reasons. In differential geometry it is the 
equation satisfied by the conformal factor z{x,y) of the metric = z^{dx^ + dy^) of a two- 
dimensional space of constant curvature [13]. In the theory of infinite-dimensional nonlinear 
integrable systems it is the prototype of a nonlinear partial differential eqnation (PDE) lin- 
earizable by a transformation of variables, involving the dependent variables (and their first 
derivatives) alone [32] 

U = 2^, = (1.3) 

In Lie theory this is probably the simplest PDE that has an infinite-dimensional Lie point 
symmetry group [35]. The symmetry algebra of the algebraic Liouville equation (1.2) is given 
by the vector fields 

X{f{x)) = f{x)dx - fx{x)udu, Y{g{y)) = g{y)dy - gy{y)udu, (1.4) 

where f{x) and g{y) are arbitrary smooth functions. 

Equation (1.4) is a standard realization of the direct product of two centerless Virasoro 
algebras and we shall denote the corresponding Lie group VIR(x) (8) VIR(y). Restricting f{x) 
and g{y) to second-order polynomials we obtain the maximal finite-dimensional subalgebra 
s[ 3;(2,M) 0s[y(2,M) and the corresponding finite-dimensional subgroup SLa;(2,M) (g) SLj^(2,M) 
of the symmetry group. 

The Liouville equation is also an excellent tool for testing numerical methods for solving 
PDE’s, since equation (1.3) provides a very large class of exact analytic solutions, obtained by 
putting 

4>ix,y) = Mx) + My): ( 1 - 5 ) 

where 4>i{x) and (j) 2 {x) are arbitrary functions on some interval I. 

In [1] Adler and Startsev presented a discrete Liouville equation that preserves the proper¬ 
ty of being linearizable and exactly solvable. In [38] Rebelo and Valiquette wrote a discrete 
Liouville equation that has the same infinite-dimensional VIR(x) (g) VIR(y) symmetry group as 
the continuous Liouville equation. The transformations are however generalized symmetries, 
rather than point ones. In onr article [23] we presented a discretization on a four-point stencil 
that preserves the maximal finite-dimensional subgroup of the VIR(x) (g) VIR(?/) group as point 
symmetries. It was also shown that it is not possible to conserve the entire infinite-dimensional 
Lie group of the Liouville equation as point symmetries. In [23] we also compared nnmerical so¬ 
lutions obtained using standard (non invariant) discretizations, the Rebelo-Valiquette invariant 
discretization [38] and our discretization with exact solutions (for 3 different specific solutions). 
It turned out that the discretization based on preserving the maximal snbgroup of point trans¬ 
formations always gave the most accurate results for the considered solutions (all of them strictly 
positive in the area of integration). 

The purpose of this article is to further explore and compare the different discretizations of 
the Liouville equation from two points of view. One is a theoretical one, namely to investigate 
the degree to which different discretizations preserve the qualitative feature of the equation: its 
exact linearizability, its infinite-dimensional Lie point symmetry algebra, the behavior of the 
zeroes of the solutions. The other point of view is that of geometric integration: what are the 
advantages and disadvantages of the different discretizations as tools for obtaining numerical 
solutions. 

In Section 2 we reproduce our previous [23] SLa,(2,M) (g) SLj^(2,M) symmetry preserving dis¬ 
cretization nsing a 4-point stencil and show that after a slight modification it can reproduce 
solutions that have horizontal or vertical lines of zeroes (or both). In Section 3 we propose an 
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alternative discretization, using a 9-point stencil, instead of the 4-point one. It approximates 
the continuous Liouville equation with precision, as opposed to the e precision of the 4-point 
discretization. We show that increasing the number of points does not allow us to preserve 
the entire infinite-dimensional symmetry algebra, nor to treat the lines of zeroes of solutions 
in a satisfactory manner. Further, in Section 3 we take a specific exact solution of the contin¬ 
uous algebraic Liouville equation (1.2) and approximate it on a 9-point lattice by a numerical 
solution. The Adler-Startsev discretization [1] is reproduced in Section 4 in a form suitable for 
numerical calculations. Section 5 is devoted to numerical tests of the invariant 4-point scheme. 
Five different exact solution of the algebraic Liouville equation are presented and then used to 
calculate boundary conditions on two lines parallel to the x and y coordinate axes, respectively. 
The solutions are then calculated numerically using four different discretizations. We compare 
the validity of the different methods and their qualitative features. Some general conclusions, 
placing this article in the context of geometric integration, are presented in the final Section 6. 

2 Point symmetries on a four point lattice 
and solutions with zeroes 

In our previous article [23] we discretized the algebraic Liouville equation (1.2) on a four point 
regular orthogonal lattice preserving the SLa;(2,M) (g) SLj^(2,M) subgroup of its Lie point sym¬ 
metries. The discretization was shown to provide good numerical results for solutions that were 
strictly positive in the entire integration region (a quadrant to the right and above a chosen 
point {xo,yo), i.e., for x > xq, y > yo). 

A particular property of the Liouville equation is that the zeroes of its solutions are not 
isolated. They occur on lines parallel to the x or y axes. Indeed, consider the infinite family 
of solutions of (1.2) parametrized by two arbitrary smooth functions of one variable </>i(x), 
4>2{y) (1-5). We take a region in which we have 4>i{x) + 4>2iy) / 0. Zeroes of u{x,y) occur 
if 4>i^x{x), or 4>2,y{y) are zero at some point Xs, or yg (or both), respectively. We then have 

u{xs,y) = 0, yy, or u{x,ys) = 0, Vx. (2.1) 

This must be reflected in any computational scheme and the value u{x, y) = 0 will also occur 
on the intersection with the corresponding coordinate axis. 

In [23] we considered several different boundary value problems. Here we restrict to the case 
of boundary conditions given on the lines x > xq, y > yo parallel to the coordinate axes. We 
can impose 

u{xs,0)=0 and/or u{0,ys) = 0 
in order to obtain a solution satisfying (2.1). 

The SLa;(2, M) (g) SLj^(2, M) invariants used in [23] to describe both the lattice and the discrete 
algebraic Liouville equation on a four point stencil were 

^ _ iXm,n-\-l Xfji,n) {Xra+l,n+l 3;m-|-l,n) 

(Xm,n Xm+l,n) (Xm,n+1 Xm+l,n+l') 


(2/m,n+l 2/m,n) 2/m+l,n) 

(2.2) 

^ ^ ; '^2 — k . 

(2.3) 

The lattice equations 


6 = 0, Vi = 0 

(2.4) 
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are satisfied by the uniform orthogonal lattice 


Xm,n = hm + xo, ym,n = kn + yo, (2.5) 

where the scale factors h and k are the same as in (2.3). The continuous limit corresponds to 
/i —>■ 0, A: —)• 0. Two further independent SLa;(2, M) ( 8 ) SLy (2, M) invariants exist on the four point 
stencil but any combination of them will either vanish, or be infinite on the lattice given by (2.4) 
(see [23]). 

The Liouville equation (1.2) was approximated in [23] by the difference scheme 
J 2 - Ji = asign(Ji)| + 6 sign(J 2 )Ji| 

+ csign(Ji)| Ji|^/^J2 + dsign(J2)| J2|^/^, ( 2 . 6 ) 

= 0 , a + b + c + d=l. 


The symbols sign(Ji) and sign(J 2 ) were omitted in [23] and were not necessary as we restricted 
our formulation to strictly positive solutions. Equation (2.6) can be solved for Um+i,n+i in terms 
of Un,m, Un+i,m and Un,m+i- On the first stencil we have m = n = 0. The boundary conditions 
are Um,o = /(nr) and uo,n = gin) with / and g given. 

Let us now rewrite the recurrence relation (2.6) in terms of Um,n, choose b = d = 0, c = 1 — a, 
a G M (in order to have an explicit scheme) and solve for Um+i,n+i- We have 


— -^m,n+l;m+l,n5 

^ 1 -|- d/l/c sign('U77^-|_i^7x't^?n,n+l) ^/lI 

■^^777. 77.~|“ 1‘7T7.~|“ 1 TL — , - 

1 + (o - l)hks\gn{Um+l,nU m,n+l) Y 


(2.7) 

( 2 . 8 ) 


The expression sign(um+i,n^tm,n+i) follows from (2.6). Here the sign before the square root is 
important since it will change when the sign of u changes in the recurrence relation. 

We shall use (2.7), (2.8) to investigate the behaviour of the numerical schemes for solutions 
that have rows (horizontal lines) or columns (vertical lines) of zeroes. We impose boundary 
conditions on the lines x = Xg and y = yg- To see the influence of the boundary conditions we 
introduce small quantities g and n on the coordinate axes that will later be set to zero. We shall 
see that these small values do not propagate elsewhere but are confined to the columns and rows 
where they were introduced. This procedure is analogous to “singularity confinement” [17, 18] 
used as an integrability criterion for difference equations. 

We first note that we have 


lim — lim ^4^ , 2 +l;m+l,n — !• 

Three cases will be considered separately: 

1. A column of zeroes. The boundary conditions are 


Umo,0 = Um,o / 0 for m ^ TUq. 

Using (2.7), (2.9) we obtain expressions for Umo,n, namely 

^mo,ri — 77- ^ 0 . 

— 1,0 

In the column to the right of the zeroes we obtain two equivalent expressions: 

'^TUQ —1,71 

Umo —1,71—1 

_ l^mo+ 1,0 

llmoH-l,n — l^mo —1,71 • 

^mo — 1,0 


(2.9) 


( 2 . 10 ) 

( 2 . 11 ) 
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Thus the zero quantity /i cancels out and Umo+i,n is finite and nonzero for all n > 0. Moreover 
Umo+i,n is expressed in terms of the given initial values and values calculated at previous nonzero 
values. 

2 . A row of zeroes can be treated completely analogously. The boundary conditions are 
replaced by 


UO,no = uo,n / 0 for n / Uq, 

and we obtain 

Um,no — l ^ 

^m,no — ?7l ^ 0 , 

"^ 0,710 —1 

i.e., a row of zeroes for = 0. The row above the zeroes satisfies 

_ ^m,no — 1 

^m,no+l — 1,710+1 5 

'^771—1,710 — 1 

_ 1 ^ 0,710 + 1 

1 ^ 771,710 + 1 — 11711,710 — 1 • 

110,710 —1 

3. Two intersecting lines of zeroes. The boundary conditions are 


( 2 . 12 ) 

(2.13) 


Wo,no = fo Umo,o = fj-; Urafl / 0 for m / mo, wo,n / 0 for n / no. 
Using the same considerations as above we find a column and a row of zeroes satisfying 

Wm,no = -fo m 7^ mo, 

WO,no —1 


Thus, for /r = 0, n = 0 the solutions Ura,n have zeroes precisely where they should. Now let us 
use (2.7), (2.9) to calculate the values of Wmo+i,n and Um^no+i-, i-e., the column at the right and 
the row above the zeroes. The final result is that (2.10) is valid for all n 7 ^ no, no + 1 and (2.12) 
for all m 7 ^ mo, mo + 1 with 

_ UmQ+l,no — l 

117710 + 1,710 + 1 — 117710 — 1,710 + 1 5 

117710 — 1,710 — 1 

while (2.11) is valid for all n 7 ^ no and (2.13) for all m 7 ^ mo. 

Finally we see that the zeroes are confined to the rows and columns determined by a zero 
in the boundary condition and that the values of Um,n everywhere else are finite, non zero and 
determined by the equations (2.7), (2.8) and the boundary conditions. In other words the rows 
and columns of zeroes do not interfere with the integration algorithm. This will be confirmed 
by numerical calculations in Section 5. 




117710 — 1,71 


7710,71 


/I 7 


n 7 ^ no, 


117710 — 1,0 
_ 117710 — 1,710 — 1 

117710,710 ~ /l^* 

110,7io — 1II 7710 — 1,0 


3 Invariant discretization of the algebraic Liouville equation 
using a larger number of points 

There are several reasons to increase the number of points on the stencil that we use. 

1. To determine whether the entire VIR(x) ( 8 > VIR(?/) symmetry group can be preserved on a 
larger lattice. 

2. To determine whether the only other SLa;(2,]R) 0 SLj^(2,M) differential invariant [23], 
namely 

I2 - g {^UUxX ^^x'} {^‘U‘Uyy 3Uy^ 


(3.1) 
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Figure 1. Points on a general lattice, e.g., xqq = x, xio = x + /iio, xqi = x + eoi; xn = x + /iio + en, 
2^20 = X + /iio + ^20 j 3^02 = X + eoi + £92, X12 = X + hiQ + eii + ei2, X21 = x + /iio + /i20 + ^21, 
X22 = x + /iio + /i20 + e2i+e22, 2/00 = y, yoi = y + koi, yio = y + Sio, yn = y + koi+Sii, 2/02 = y + koi + ko2, 
2/20 = 2 / + <^10 + <^201 2/12 = 2 / + ^01 + ^02 + ^12, 2/21 = 2 / + ^01 + ^10 + <^20) 2/22 = 2 / + ^01 + fco2 + 5i2 + 622- 


l0.2) . 

b.2) 

[2,2) 

Ill 

IV 

[2,^) 

1 

i'm 1 

II 

.1:2)- , 

^2.0) 


Figure 2. A stencil for the 9-points scheme. 

can be invariantly discretized on a larger lattice. Four points are clearly not sufficient to 
approximate two first- and three second-order derivatives. 

3. To approximate the algebraic Liouville equation with a higher degree of accuracy in h 
and k and thns possibly improve the numerical calculations. 

In Section 2 and in [23] we have shown that the Liouville equation can be approximated on 
4 points. To approximate an arbitrary second-order PDE for a fnnction u{x,y) we need at least 
6 points. An invariant discretization may need more than six. 

Equation (2.6) satisfies 

^hrn^ ^ { J 2 - Ji - [a| Ji |3/2 + bJi\ J 2 \+ c| Ji|J 2 + d\ J 2 } 

— UxUy U j [1 “h 

and thns provides a first-order approximation of the algebraic Liouville equation. In this section 
we will explore a second-order approximation (order 0{h‘^,k'^,hk)) of the equation (1.2). To 
do this we shall use a 9-point stencil as shown on Figs. 1 and 2. The 4 well-behaved inva¬ 
riants (2.2), (2.3) make nse of the four vertices of rectangle I on Fig. 2. Instead of the vertices of 
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rectangle I we could use any other 4 points, and we shall use the vertices of the rectangles II, III 
and IV. The invariants involving the independent variables and r]a (a = 1,... ,4) all vanish 
on the orthogonal lattice (2.5). The invariants depending on the dependent variables Uij that 
are finite and nonzero on this lattice are 


Ji = uoiuioh^ , 

Js = uiiU2oh‘^k‘^, 
J 5 = uiiUo2k?k^i 
Jj = Ui2U2lh?‘k^, 


J 2 = uoouiih^ k‘^ , 
J4 = uioU2ih‘^k‘^, 
Jq = UQiUi2h‘^k‘^, 
Jg = UiiU22h^k‘^- 


The quantities Ji,..., Jg are linearly independent but one polynomial relation exists between 
them, namely 


J4 Je 


Ji J 7 . 


The continuous limit is obtained by expanding the invariants into Taylor series and then taking 
—)■ 0, A: —)• 0. We shall assume that they tend to zero at the same rate, i.e., k = ah, a ~ 1. 
We have 

Ji = + huux + kuuy + hkuxUy + ^h'^uuxx + \k‘^uuyy + •••], 

J 2 = h?k^ + huUx + kuUy + hkuuxy + \h?‘uuxx + \^‘^uuyy + •••], 

T 3 — h k T ShuUx T kuUy -\- 2 k hki^uuxy T ‘^rLyVLx') ~\~ h {^nuxx T T * * ■ j, 

T 4 — h k j^zx T ShuUx T kuUy ~\~ 2 k ^^yy T hk^^unxy T UyUx^ ~\~ h {^nuxx T T * * ■ j, 

J 5 = h^k^ [zt^ + huUx + "ikuuy + k^{^uuyy + 2 zty) + hk{uuxy + 2uyUx) + \h?‘uuxx + •••], 

Je = h?k^ [u^ + huUx + "ikuuy + k^{^uuyy + 2zty) + hk{ 2 uuxy + UyUx) + \h?uuxx + •••], 

J 7 = h'^k'^ [u^ + 2>huux + Skuuy + k‘^{2Uy + ^uuyy) + hk{^uuxy + bUyUx) 

+ ( 2 zz^ + ^uuxx) +■■■]) 

Jg = [zz^ + Shuux + 3 kuuy + k'^{^uuyy + 2zty) + hk{ 5 uuxy + AuyUx) 

+ (| ziZZa ; a ; + 2 zz ^) + • • • ]. ( 3 - 2 ) 

We see that U 22 , U 02 , U 20 and zzqo figure only once each in the invariants, namely in Jg, J 5 , J 3 
and J 2 , respectively. On the other hand zzqi, uiq, U 12 , and U 21 figure twice each, respectively in 
(Ji, Jq), (Ji, J 4 ), (Jg, J 7 ) and (J 4 , J 7 ). The value zxn figures in all four of J 2 , J 3 , J 5 and Jg. 

To lowest-order we have 


J2 Ji — J4 J3 — '^6 '^5 — '^8 J 7 — h k (^uUxy UxUy') [I T 0(^h, A^)]. 


To obtain the left hand side of the algebraic Liouville equation (1.2) up to order Oih?, hk, k'^) 
we need the differences J 2 a — J 2 a-i to a higher-order than in (3.2), namely 


J 2 Jl — h k •( {uUxy UxUy^ T {nUxxy '^y'^xx^ T ^ i^'^'^xyy '^x'^yy') ^ ? 


h 


k 


^xy ’ 2 ujujxxy ^y^xx) ' ^y^^xyy 


Ja Js — h h UxUy') ^ {^^xxy UyUxx^ H“ ^{^"^xyy ? 

- Js = hV - u^u,) + ljuu.,, - - «.«„)} . 

e/g Jj — h h “1“ ~^{.^^xxy ^y'^xx) ~^{.^^xyy • 


(3.3) 
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In [23] equation (1.2) was approximated to order 0{h,k). To approximate it to 0{h‘^,hk,k‘^) 
we must get rid of the terms of order 0{h, k) in (3.3). 

The left hand side is approximated to the needed order by 


«[4(J2 - Ji) - (Je - J 5 + J4 - h)] + mh - J 7 ) - 3(J6 -J 5 + J 4 - J 3 )] 
= 2h^k^{uuxy — UxUy){a — /3) [l + 0[h^, hk^ /c^)], 


(3.4) 


where a and /? are arbitrary real constants. 

To express the right-hand side of (1.2) we use the basis 

Bi = — Jg) = h^k^vk^{\ + i?i), B2 = J4 + Jq — Jg — Ji = h?k‘^R2, 

Bg = J2 — Ji — h^k'^Rg, B4 = J4 — Jg — k'^ R4, Bg = Jq — Jg — h^k'^Rg, 

Bq = Jg — T7 = h'^k'^Rg, ( 3 . 5 ) 


where Ri,... ,Rg are all of the order 0{h'^, hk, k‘^). 

The left hand side of (1.2) is already expressed in this basis (see (3.4) using Bg,... ,Bq). 
From the basis elements (3.5) we can calculate v? as 


Bi -|- ^ CiBi = h'^k'^v? [1 + 0{h^, hk, /c^)]. 


i=2 


with 5 free real parameters q. To obtain we have several possibilities. One is to take 
/ 6 \ 3/2 

^ CiBi = h^k^u^ [1 + 0{h‘^, hk, k"^)]. 


(3.6) 


The corresponding discrete Liouville equation is then 


a 


[4^3 - (^4 + Bg)] + /3[4S6 - 3(^4 + B 5 )] = ^ CiBi 


i=2 


Bi + CiBi 


i=2 


1/2 


(3.7) 


with 


2(a-/3) = 1. 


(3.8) 


Another possibility is to replace the basis (3.5) by 
Ai = Bi, Aa = Bi + -^Ba, a = 2,... ,6. 

We can then approximate the right-hand side of the discrete algebraic Liouville equation by 

6 6 

7a,bAaV\M = la,b[l + 0{h\hk,k^)]. 

a,6=1 a,6=1 

Then the discrete algebraic Liouville equation reads 


6 

2a[4A3 - 2Ai - (A 4 + A 5 )] + 2/3[AAg + 2Ai - 3(^4 + Ag)] = ^ ia,bAa^/\M, (3.9) 

a,6=1 

6 

la,b = 2(a - 13). 

a,6=1 
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The general invariant equations (3.9), (3.6) use all 9 points on the stencil and contain a lot of 
free parameters. The parameters can be chosen to simplify calculations, though the choice of 
(a, /3, Cl,..., cg) or (a, /3, ^a,b) is restricted by the type of boundary conditions we wish to impose. 

The quantity U 22 figures in Jg only. An explicit scheme is obtained if Jg figures linearly in 
the corresponding invariant discrete Liouville equation. 

One possibility is to choose a = —3/3 and 02 = 03 = 04 = 05 = 0, cg = 1/2 in (3.7), (3.8). 
Then /3 = —1/8 and the invariant Liouville equation reduces to 




J 7 + 3 (J 2 — Ji) 


^(3Ji — — J^\. 


(3.10) 


In terms of the field Uij (3.10) reads 


U22 — - [U12U21 + 3 (uiit(oo — UiqUoi) 

un 


1 

71 


hk{3uoiUio - Ui2U2l)\/\3uoiUio - Ui2U2l\]. 


(3.11) 


so that U 22 is expressed in terms of uqq, uqi, uiq, uh, uu and M 2 I) i-e., only 7 points are involved. 

Another simple possibility is to choose a = —3/3 and 7^5 = in (3.9). Then we have 

[3 = —1/8 and we obtain 


J 7 + 3(J2 — Ji) — sign(3Ji — J7)^y\3Jl — Jy\ 
1 - ^ sign(3Ji - J 7 ) v^|3Ji - J 7 I 

In terms of the field Uij (3.12) reads 


(3.12) 


U22 = \ U12U21 + 3 {unUoo - uiouoi) 


- -^hkuoiUiosign{3uoiUio - ui2U2i)y/\3uoiUio - ui2U2i\ 


X < Uii 


hk 


1 - ^ sign(3noittio - ui2U2i)y/\3uoiUio - ui 2 U 2 i\ 


-1 


(3.13) 


Again only 7 of the 9 points on a stencil are used. 

Equations (3.11) and (3.13) are to be viewed as recursion relations, expressing M 2,2 in terms 
of 6 points on a rectangle of which the point (2,2) is the top right vertex (see Fig. 2). 

By construction (3.11) and (3.13) are better approximations of the equation (1.1) than is (2.6). 
This does not mean that they will provide better numerical results and some comments are in 
order. 

1. Boundary conditions for a numerical solution on a 4-point lattice require the knowledge 
of u{x,y) on two lines, e.g., Umfl and Mo,nj i-e., m(x,0) and u{0,y). On the 9-point lattice we 
must start with 2 sets of parallel lines, e.g., Umfi, Um,i and Mo,n, ui,n- This amounts to giving 
m(x, 0), M(0,y) and the first term of Uy{x,0), Ux{0,y)- This is more information than is needed 
in standard (non invariant) discretizations and indeed more information than is needed in theory 
to determine a solution completely. Hence once u{x, 0) and m(0, y) are given Uy{x, 0) and Ma;(0, y) 
cannot be chosen arbitrarily. In our numerical solutions we calculated the conditions using an 
exact solution so m(x,0), u{0,y) and Ux{0,y), Uy{x,0) are consistent. 

2. Contrary to the case of a 4-point lattice, instabilities close to zero lines of solutions cannot 
be avoided on 7- or 9-point lattices. Indeed let us give initial conditions on the first square 
satisfying mqo / 0, mqi / 0, miq / 0, mh / 0, M12 / 0, U 20 = ei, M21 = € 2 - From the known 
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solution of the PDE (1.1) we expect the solntion to satisfy U 2 ,m = 0 for m >2. Equation (3.11) 
implies 

1 

U22 = - 

Mil 

Thus U 22 is not strictly zero for ei = e 2 = 0, it does however satisfy U 22 rsj 0{h‘^,k^,hk). This 
is acceptable, however the problem arises when we shift the stencil and calculate M 32 which is 
supposed to be finite and nonzero if we assume M 30 / 0, M 31 / 0. What we obtain from (3.11) is 


Ui2e2 + 3(miiM 00 - ^^10^i0l) - -^hk{3uoiUio - Mi2e2)\/|3MoiMio - Mi2e2| 


- -^hk(3eiuii — 
v2 


Thus, M 32 is singular for 62 = 0 and becomes finite only in the continuons limit h = k = 0. This 
will quite obviously create numerical instabilities. They are avoided only for very special initial 
conditions, such that U 22 = 0 for all h and k. Using (3.13) leads to the same kind of problems. 

Sadly (for the 9 points scheme) the answers to all three questions posed in the beginning of 
this section are negative. 

1. The only VIR(x) ® VIR(?/) invariant lattice is given by requiring conditions (2.2) in all 
fonr quadrangles of Eig. 2. Indeed we have 


prX(x^)Ca = (xii - xoo)(a:io - 

pvY{y^)r]a = {yii - yoo){yio - yoi)ila, a = 1 , 2 ,3, 4 


so (weak) invariance requires = Va = 0 and adding further points does not help. Moreover, no 
function of Ji,..., Jg has the correct continuous limit and is invariant under VIR(x) ( 8 > VIR(y). 

2. We have determined that the invariant (3.1) cannot be discretized in an SLa;(2,]R) ( 8 ) 
SLy(2,M) invariant manner nsing Ji,..., Jg but we do not present the details here since this 
question is not related to the Lionville equation. 

3. Nnmerical results for several exact solutions showed that serious instabilities occur for the 
7-point scheme. As one can see in Eig. 3 representing the solution fi given in (5.9), instabilities in 
the 7-point case occur almost immediately. The same is true for other solutions of equation (1.2). 

Our conclusion is that the 9-point (or 7-point) scheme is too unstable to be useful. We present 
it here because we think that this negative result is not a priori obvious and that this discussion 
may be usefnl. 


4 The Adler—Startsev linearizable discrete Liouville equation 


Adler and Startsev [1] have presented a discretization of the algebraic Liouville equation (1.2) 
on a four-point lattice, namely 


®m+l,n-|-l (It 


1 


1 + 


1 


= 1 . 


^m+l,ri. / \ ^m,n+l, 

This equation is linearizable by the substitution 

^m,n) 


(4.1) 


Ojm.r). — 


where bm,n satisfies the linear equation 

^m+l,n+l T — 0* 
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Figure 3. Comparison between the 4-point and 7-point approximation for the solution /i given in (5.9). 
The dark grey graph represents the results of the application of (3.12), (3.13), while the light grey one is 
obtained by using the 4-point formula (2.7), (2.8). The integration is carried out starting from the the 
point (—1.0, —1.0), the recursion formula is applied for steps h = k = 0.02 on a grid 20 x 20. 


Hence the exact general solution of (4.1) is 

_ {Cm+l Cfn){kn+1 ^n) 

{Cm+1 + kn){Cm + kn+l) 

where Cm, k^ are arbitrary functions of one index each. 

In [23] we showed, following [1] that the continuous limit of (4.1), for am,n = —^Um,n 
when h and k go to zero, gives (1.2) and that it has no continuous point symmetries but must 
have generalized symmetries. Moreover by defining Cm = kn = (p 2 {ym,n) with Xm,n 

and ym,n defined in (2.5) we have 

Cm+l = (t)l{x) + -b 0(/l^), kn+l = 4’2{y) + k^^ +0{k‘^), 

and thus am,n = + 0{h^, h‘^k, hk"^, k^) a first-order approximation of the general 

solution of (1.2) given by (1.3) . 


5 Numerical tests of the 4-point scheme 


In this section we shall apply the invariant recursion formula (2.7), (2.8) to solve a set of 
boundary value problems on a quadrant in the xy-plane. Boundary conditions will be given on 
two orthogonal lines parallel to the x and y axes, respectively, and numerical solutions will be 
constructed above and to the right of these lines. The numerical solutions will be compared with 
exact solutions of the continuous equation for the same boundary conditions. In practice we will 
start from exact solutions given by choosing 4>{x,y) = 4>i{x) + </> 2 (y) in (1.3) and calculate the 
values of these functions on the boundaries. The global estimator which we use is the discrete 
analog of relative distance in L|,. We compute the quantity 

where Fij are the values of the exact solution F on the lattice sites and T)", with a = Inv, AS, RV, 
or stand are the values computed numerically for the invariant, Adler-Startsev, Rebelo-Vali- 
quette or standard discretization, respectively. A similar analysis is performed for the other 



(5.1) 
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recursion formulae. The summation will be over all points of the lattice for which the calcula¬ 
tion was performed. 

The Xa quantity (5.1) provides information about the overall averaged behavior of numerical 
solutions, rather then about their point-by-point behavior. Geometrical features of solutions are 
better reflected by plots of individual solutions and by a relative error function such as 


Rij — 


— F- • 


F, 




(5.2) 


A characteristic property of all solutions of the algebraic Liouville equation concern the zeroes. 
They either have no zero in any finite domain F, or the zeroes are not isolated, but occur along 
continuous lines parallel to the x or y axes. 

We will compare results using four different discretization methods and thus four different 
recursion formulae, expressing Um+i,n+i hr terms of Um,n, Um+i,n and Um,n+i- For comparison 
we present the four formulae for the first position of the stencil, i.e., m = n = 0. In all cases 
the left hand side of (1.2) is approximated by 


unuoo - uiouoi = h^k^{uuxy - u^Uy), 


where h and k are the lengths of the steps in the x and y directions, respectively. The right-hand 
side of (1.2) is approximated differently in each case. The corresponding recursion formulae and 
their continuous limits up to one order beyond the leading one are: 

1. The invariant method (2.7), (2.8) (preserving the SL2,(2,]R) (g) SLy(2,]R) symmetry group 
as point symmetries) 


unuoo - uiouoi = hk[auoiuio + (1 - a)uooMii] sign(uoiitio)\/|'«oi«io|, 


(5.3) 


k h 

UUxy UxUy — U -\- (jiu Uy Uyyllx UUxyy^ ^ “t” (jiu Ux “t” UyV>xx UUxxy^ (^•^) 

2. The Rebelo and Valiquette method (preserving the entire infinite-dimensional symmetry 
algebra as generalized symmetries) 


unUoo — — hkuoouoiuio, 

UUxy UxUy — U -\- (^‘2u Uy UxUyy UUxyy^ ^ {^U Ux “t” UxxUy UUxxy^ 

3. The Adler-Startsev method (preserving linearizability of the Liouville equation) 

^toi + Uio 


— uiouoi — hkuooun 


— hkuoiuio 


UUxy UxUy — U “b (2'U Uy -\- UxUyy UUxyy^ ^ T ^2u Ux “b UxxUy UUxxy') 2 ■ 


(5.5) 


(5.6) 

(5.7) 


4. The standard method (not preserving any specific strncture) defined on the 4 points of 
a square lattice is 


unUoo - uiouoi = kku^Q, (5.8) 

o . . h . k 

UUxy UxUy U “b [UxxUy UUxxy) “b [UxUyy UUxyy) ~‘ 

Each of these formulae gives a different explicit expression for un in terms of the already 
known values of uoo, and uio- 
Several comments are in order: 
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Table 1. Relative mean square distance (5.1) between the numerical solutions and the analytic one in 
the domain Vq. 



Xinv 

XAS 

Xrv 

Xstand 

fl 

5.2 X 10“® 

2.7 X 10-5 

3.1 X 10-^ 

9.2 X 10-1 

f2 

3.4 X 10"^ 

1.5 X 10-^ 

7.6 X 10-3 

2.2 X 10-2 

h 

4.7 X 10-® 

1.5 X 10-5 

3.0 X 10-3 

9.2 X 10-3 

h 

4.3 X 10-5 

7.9 X 10-5 

5.2 X 10-3 

2.0 X 10-2 

h 

3.8 X 10-2 

3.0 X 10-2 

2.8 X 10-1 

4.3 X 10-1 


1. In (5.4) there is no dependence on the parameter a. It will only appear at the order hk 
(not /i^ or k'^). That is the reason why the dependence of the numerical results depends weakly 
on the choice of a (see Fig. 7 below to confirm this). 

2. Formulas (5.4) (invariant method) and (5.7) (Adler-Startsev method) coincide. They differ 
in the higher-order terms. Table 1, as could be expected, confirms that the results obtained by 
these two methods are similar. 

We consider 5 different solutions of the continuous algebraic Liouville equation (1.2), namely 


/i 

/2 

/3 

U 

h 


2 __ 

(x^ -h 1) (y2 -h 1) (tan“^(x) -|- tan“^(y) -|- 6)^ ’ 

8 (l - 4 (x -h 5)) (1 - 4y) exp (-2x(l -h 2x) - 2y{y -h 2)) 

(^^-2x(l+2x) g2j/(l-2y) ^^2 

3.38 sin(1.3(x -|- 10“^)) cos(1.3(y -|- 10“^)) 

(cos(1.3(x -h 10“2)) + sin(1.3(y -|- 10“2)) + 3)2 ’ 

8xy 

(x2 -h y2 _p 2 )^’ 

383 le3-862(2.5(a:-0.5)-|-0.4y-|-2.5) 

('g9.655(a:-h0.5) + 12.8‘ie^-^^^y)^' 


(5.9) 


The functions fi and /s do not contain any zeroes in any finite domain. The functions /2 
and /4 have one row and one column of zeroes each. The function /s contains infinitely many 
orthogonal lines of zeroes, since it is a periodic function. Finally, /s is a wall like function, with 
no zeroes. We mention that for /i, /s and /4 the first-order corrections in (5.4) and (5.7) vanish. 
Plots of the exact solutions /i, • • ■, /4 are given in Fig. 4, /s on Fig. 6a below. The right-hand 
sides of (5.5), (5.6) and (5.8) are polynomials whereas the invariant case (5.3) involves square 
roots. 

The numerical computations were performed on the square domain Dq = [—1.5,1.!] x 
[—1.0,1.6], with steps of equal length h = k = 0.02, for a lattice of 130 x 130 points. Somewhat 
arbitrarily we choose the parameter a in the symmetry invariant recursion formula to be a = 1.0. 
The boundary conditions are given on the bottom and left side of the square. 

In Table I we give the y quantity (5.1) for all five solutions using 4 different methods. We 
see from Table 1 that xinv and xas are in general of the same order, as are XRV and Xstand- 
The values of xinv and XAS are better than those of XRV and Xstand by at least one order of 
magnitude, usually by 2 orders, with xrv always better than Xstand- The faster the solution 
changes, the worse is the result (for all methods), specially for the solutions /2 and /s. 

In order to test the stability of the algorithms with respect to the size of the adopted meshes, 
we made another series of calculations involving the above test functions over a fixed domain 
Vi = [—1.905,1.895] X [—1.905,1.895], larger than Vq, and spanned it using different lattice 
scales with h = k. A general flavor of such calculations can be extracted from Table 2 for 
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Figure 4. Plots of the exact solutions (a) /i, (b) / 2 , (c) /a and (d) fi in the domain Vi used in the 
numerical integrations. 


Table 2. Log^pX in the domain Vi for each discretization procedure we considered changing the size of 
the lattice meshes. We report only the results for /i with the parameter a = 1. 


h = k 

Xinv 

XAS 

XRV 

Xstand 

8.X 10“^ 

-3.98219 

-4.30917 

-2.57467 

-2.13488 

4. X 10-1 

-4.5782 

-4.91072 

-2.87464 

-2.41647 

2. X 10-1 

-5.17707 

-5.51254 

-3.17517 

-2.7076 

1. X 10-1 

-5.77759 

-6.11449 

-3.47595 

-3.00363 


the function /i For the solution fi (with no zeroes) the value of xinv and XAS are comparable 
and at least two orders of magnitude lower than the other two. The values of xrv are always 
lower than Xstand but of the same order. Generally speaking, decreasing the scale of the mesh 
by a factor of 0.3 implies decreasing the value of x by a factor of 0.6 for the invariant and 
Adler-Startsev discretization and by a factor of 0.3 for the other two discretization (qnadratic 
as opposed to linear convergence). 

Let us now turn to the point-by-point behavior of the solutions. As discnssed in Section 2 
a characteristic property of all solutions of the Lionville eqnation is that zeroes are not isolated 
but occur in straight lines parallel to the axes. To see how well this is reflected in numerical 
solutions, let us first concentrate on solution /2 which has zeroes on horizontal and vertical 
lines passing through the saddle point x = —0.25, y = -1-0.25. These points are not on the 
lattice due to the definition of the domain Vi with h = k = 0.01. On Table 3 we give the 
values of the solution /2 at the four lattice points nearest to the saddle point. The AS solution 
has the first 4 digits coinciding with the exact one, the invariant one has 3. The RV and 
the standard method have 2. Similar results are also valid for the other solutions with zeroes 
ih and /4). 
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a) 




b) 


d) 



Figure 5. Plots of the relative error Rij defined in (5.2) for the numerical approximation of the solution /2 
in the domain Vi using the individual discretizations: (a) Invariant, (b) Adler-Startsev, (c) Rebelo- 
Valiquette and (d) standard. The graphs are not in the same scale. The maximal value of the error is 
approximately 2 x 10“^, 2 x 10“"*^, 1.5 x 10“^, 4 x 10“^ for the discretizations (a), (b), (c ) and (d), 
respectively. 


Table 3. Value of /2 at the four lattice points nearest to the saddle point (—0.25, +0.25) for the exact 
solution and the four numerical approximations. 


point 

Exact 

Inv 

AS 

RV 

stand 

(-0.255,0.255) 

(-0.245,0.255) 

(-0.255,0.245) 

(-0.245,0.245) 

-4.14391 • 10“"“ 
4.14391 ■ 10-* 
4.14391 ■ 10-* 
-4.14391 • 10'^ 

-4.1419 • lO""' 
4.14257- 10"^ 
4.14123 • 10"^ 
-4.1419 • lO"'^ 

-4.14369 • lO"'^ 
4.14369 • lO”'^ 
4.14369 • 10“^ 
-4.14369 • 10"'‘ 

-4.1281 • lO"'^ 
4.12877- 10-* 
4.12744 - 10-* 
-4.1281 - lO""^ 

-4.11579 - 10"^ 
4.11513 - 10“'‘ 
4.11645 - 10-^ 
-4.11579 - 10-^ 


The behavior of the solution /2 is plotted on Fig. 5 for the entire region Di where we show 
the values of the error function Rij of (5.2). The maximal value of Rij is 2 x 10“^ in Fig. 5b 
for the AS method, 10 times higher for the invariant method, 75 higher for RV and 200 times 
higher for the standard one. For the invariant method the error is concentrated at the saddle 
point (see Fig. 5a) with a tail behind it. The AS method has maximal error at the four extrema 
(see Fig. 5b) with no tail. The other two methods have maximal errors on the maxima of the 
solutions (not the minima) with tails in both directions. 

On Fig. 6 we analyze the solution /s in detail in a point-by-point manner in the domain Ri. 
Fig. 6a represents the exact solution (a wall of constant height). It has no zeroes anywhere in the 
finite real plane. The height of the wall gradually decreases for solutions b, d and e, but much 
more slowly for the invariant method b. For the AS method the height increases so we present 
the height of the solution on Fig. 6c in a different scale. The increase in Fig. 6c is comparable 
with the decrease in Fig. 6b (a factor of about 2.5). 

Finally we analyze the role of the parameter a in the formula (2.7), (2.8). The continuous 
limit (5.4) shows that a appears for the first time in terms of order 0{hk). Numerical calculations 
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Figure 6. Pictures of the function in the domain Vi for the exact solution (a) and for the individual 
discretizations: (b) Invariant, (c) Adler-Startsev, (d) Rebelo-Valiquette and (e) standard. It is clear by 
inspection that in this case all proposed discretizations are numerically unstable, leading to decaying or 
blow up of the computed solutions by recursive formulae. 


of xinv for the function /2 shows a variation of about a factor 2, when a G [—0.5,1.5] (see Fig. 7). 
Furthermore, this function takes a minimum for a = 0.17. However, this value is strongly 
dependent on the test function considered. 


6 Conclusions 

Both from the point of physics and from the point of view of geometric integration we see that for 
discretizing the Liouville equation we have to choose which characteristic feature of the equation 
we wish to preserve. Adler and Startsev [1] have shown how to preserve linearizability and the 
existence of a class of exact solutions depending on two arbitrary functions of one variable. We 
have shown that for a wide class of solutions a recurrence formula based on their method provides 
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Figure 7. Plot of xinv for the function /2 in the parameter range a G [—0.5,1.5]. 


the most accurate results, both using the global criterion and comparing local point by point 
convergence using the Rij criterion. On the other hand, linearizability, just like integrability is 
a property of a very restricted class of nonlinear PDEs. 

The existence of a nontrivial Lie point symmetry group is a much more generic property, 
specially for PDEs coming from fundamental physical theories. Erom this point of view the 
Liouville equation is again special: its Lie point symmetry group is infinite-dimensional. Rebelo 
and Valiquette [38] have presented a discretization that preserves this entire infinite-dimensional 
symmetry group as a special type of generalized symmetries. As opposed to more general higher 
symmetries, their symmetries have a global group action and are very interesting from the 
theoretical point of view. Erom the numerical point of view of we have shown that that the 
precision of the RV solutions is systematically better than that of those obtained by the standard 
method (though of the same order of magnitude). The measure of the validity is given by the 
quantity x of (5.1). 

Einally, the method proposed in [23] and further developed in this article preserves point 
invariance under the maximal finite subgroup of the infinite-dimensional symmetry group. Nu¬ 
merical methods based on this partial preservation of symmetries perform very well for all 
solutions in some case even better than the Adler-Startsev case. 

In future work we plan to study symmetry preserving discretizations of other equations with 
infinite-dimensional symmetry groups, such as the Kadomtsev-Petviashvili equation, and the 
three-wave interaction equation. A symmetry preserving discretization of the Korteweg-de Vries 
equation has provided encouraging results [3]. 
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